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ABSTRACT 

We present a detailed study of the internal kinematics of the Galactic Globular 
Cluster M4 (NGC 6121), by deriving the radial velocities from 7250 spectra for 2771 
stars distributed from the upper part of the Red Giant Branch down to the Main 
Sequence. We describe new approaches to determine the wavelength solution from day¬ 
time calibrations and to determine the radial velocity drifts that can occur between 
calibration and science observations when observing with the GIRAFFE spectrograph 
at VLT. 

Two techniques to determine the radial velocity are compared, after a qualitative 
description of their advantages with respect to other commonly used algorithm, and 
a new approach to remove the sky contribution from the spectra obtained with fibre- 
fed spectrograph and further improve the radial velocity precision is presented. The 
average radial velocity of the cluster is (v) = 71.08 ± 0.08 kms“^ with an average 
dispersion of = 3.97 kms“^. Using the same dataset and the same statistical 
approach of previous analyses, 20 additional binary candidates are found, for a total 
of 87 candidates. A new determination of the internal radial velocity dispersion as a 
function of cluster distance is presented, resulting in a dispersion of 4.5 kms“^ within 
2' from the center of cluster and steadily decreasing outward. We statistically confirm 
the small amplitude of the cluster rotation, as suggested in the past by several authors. 
This new analysis represents a significant improvement with respect to previous results 
in literature and provides a fundamental observational input for the modeling of the 
cluster dynamics. 

Key words: globular clusters: individual: NGG 6121 - stars: kinematics and dynam¬ 
ics - techniques: spectroscopic - techniques: radial velocities. 


1 INTRODUCTION 

M4 (NGC 6121) is one of the most studied Galactic Glob¬ 
ular Clusters (GGC). Its proximity, brightness and the 
moderate concentration make possible the study of its 
inner regions in great detail. This cluster has been the 
object of several intensive photometric campaigns, both 


* Based on observations collected at ESO Paranal Observatory 
within the observing program 71.D-0205, 77.D-0182 and 383.D- 
0802 
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from space (e. g. Bedin et al. 2013, Paper I, Piotto et al. 
2015) and from the ground (e. g. Kaluzny et al. 2013, 
Libralato et al. 2014). Red giant members of this cluster 
are accessible to high-resolution spectroscopy with even 
medium-size telescopes, resulting in a wealth of abundance 
studies (e. g. Brown & Wallerstein 1992, Ivans et al. 1999, 
Smith & Briley 2005, Yong et al. 2008, Marino et al. 2008, 
D’Orazi et al. 2013), which have also provided evidences of 
multiple populations. Only recently the presence of multiple 
populations has been photometrically observed in the Red- 
Giant Branch (RGB) (Monelli et al. 2013) and in the Main 
Sequence (MS) (Milone et al. 2014, Paper II, Nardiello et al. 
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2015) This cluster also represents a natural benchmark for 
many theoretical studies, such as numerical simulations of 
the internal dynamics (Heggie & Giersz 2008; Heggie 2014) 
and the interactions with the Galaxy (e. g. Dinescu et al. 
1999). 

Despite being so well-studied, the internal dynamics of 
this cluster are still poorly known. The kinematics of other 
GGG’s have been explored in more detail. One purpose of 
such studies has been to determine the presence of a cen¬ 
tral intermediate-mass black hole in a cluster center, us¬ 
ing integral-held spectroscopy to measure the radial velocity 
dispersion in the center of the cluster (e.g. within 10”; see 
for example Feldmeier et al. 2013, Liitzgendorf et al. 2013), 
Theories of gravity have also been tested by measuring ra¬ 
dial velocities of individual stars in the outer regions of a 
cluster (farther than 1' from the center of the cluster, see 
for example Baumgardt et al. 2009, Ibata et al. 2011). Both 
techniques have their disadvantages when determining the 
kinematical status of a cluster: integral held spectroscopy 
requires a deconvolution of the observed spectra along the 
line of sight, while results from individual stars may be ham¬ 
pered by the low statistics, especially near the center of far 
clusters where it may hard to observe individual stars. 

M4 represents an excellent target for kinematic analy¬ 
sis, thanks to its proximity and the many detailed studies of 
its properties. Still, the latest determination of the internal 
radial velocity dispersion as a function of distance dates back 
to the radial velocity study of Peterson, Rees & Gudworth 
(1995, hereafter P95) , which sampled 182 members with 
radial velocity precision 1 kms“^. Despite the many recent 
spectroscopic campaigns on M4, more precise and accurate 
wavelength calibrations are required for radial velocity de¬ 
termination than for chemical analyses. Spectra obtained 
for abundance studies may be ahected by systematic shifts 
in radial velocity among different observations, or simply 
lack the required spectral resolution for a precise determi¬ 
nation of velocities. These aspects are usually not important 
for chemical analyses, but they severely impact the use of 
spectroscopic data from different instruments for detailed 
kinematics analyses. 

During last decade our group has collected thousands 
of spectra using the GIRAFFE spectrograph at VLT to es¬ 
timate the geometrical distance, by comparing the radial 
velocity dispersion with the proper motion dispersion, and 
dynamical state of several GGC. Analyses of these spectra 
have been hampered by systematics in data calibration and 
radial velocity measurements; see for example Milone et al. 
(2006) and Sommariva et al. (2009; hereafter S09). 

This project is the spectroscopic counter-part of an ex¬ 
tensive astrometric and photometric investigation of M4, 
involving about 700 WFG3/UVIS and 120 AGS/WFG HST 
images (GO-12911, PI Bedin), aimed at searching for and 
characterizing the binary population in the core as well as 
in the outskirts of the cluster (Bedin et al. 2013, Paper I for 
a careful description of the HST large program GO-12911). 

In this paper we have focused our attention on signifi¬ 
cantly increasing wavelength calibration accuracies and ra¬ 
dial velocity determinations of very low SNR stellar spec¬ 
tra. After introducing the dataset in §2, we devote §3 to a 
detailed description of the improvements in wavelength cal¬ 
ibration. Our methodology to properly account for the ge¬ 
ometrical distortions in a CGD chip when determining the 


drift in radial velocity that may occur between day-time 
calibration and science observations is described in §4. In §5 
three well-known techniques are compared in order to estab¬ 
lish which one is the most suitable to the characteristics of 
the instrument. 

Finally, in §6 we discuss the application of these im¬ 
proved data reduction techniques to a dataset comprising 
7250 spectra obtained with GIRAFFE for 2771 stars of the 
GGG M4 in the wavelength region 5145 — 5360 A. Only 
about 300 stars in our sample reside on the RGB. The vast 
majority of the stars are members of the Sub-Giant Branch 
(SGB) and MS stars, i. e. relatively faint targets with low 
SNR spectra. Careful consideration of the continuum place¬ 
ment and the sky contribution for each stellar spectrum is 
required. These aspects have already been discussed in Mala¬ 
volta et al. (2014; hereafter LM14). While in our previous 
work we focused on the determination of atmospheric pa¬ 
rameters over a wide range of evolutionary states and spec¬ 
tra of different quality, in this paper we aim to increase the 
precision of the radial velocities available in the literature 
and provide new insights into its internal kinematics. 


2 THE SPECTROSCOPIC DATA SET 

The spectra for our project were originally used by S09 to 
study the internal velocity dispersion of M 4 and to search 
for spectroscopic binaries. As described in LM14, the instru¬ 
mental configuration was the GIRAFFE medium-high reso¬ 
lution spectrograph fed by the VLT Fibre Large Array Multi 
Element Spectrograph (FLAMES; Pasquini et al. 2000) in 
MEDUSA multi-hbre mode. This configuration produced 
single-order spectra for 132 objects (target stars and sky) 
in each integration. Using the GIRAFFE HR9 setup yielded 
spectra with dispersion of 0.05 A/pixel and measured 4-pixel 
resolving power R = A/AA = 25,800 in the wavelength range 
5143 A < A < 5356A. The HR9 wavelength region is partic¬ 
ularly useful for radial velocity studies because it contains a 
large number of iron lines and very strong absorption lines, 
notably the Mg-b triplet. 

The M4 target stars were selected by S09 from an 
astrometric and photometric catalog based on Wide Field 
Imager (WFI) data from the ESO/MPIA 2.2m telescope 
Anderson et al. (2006). Each fibre’s radius on the sky is 0.6”, 
so each chosen star was constrained to have no bright neigh¬ 
bors (defined as Vneighbor < Vtarget + 2.5) within an angular 
distance of 1.2”. 

The spatial distribution and photometric properties of 
our target stars are the same as illustrated in Fig. 1 and 2 
of LM14. Our target stars essentially sample all of M 4 spec¬ 
troscopically up to twice the half-mass radius. The limiting 
magnitude of the targets, V < 17.5, was chosen to ensure 
that a single M4 integration had signal-to-noise SNR > 10 
for each star, which was considered by S09 to be the mini¬ 
mum threshold to determine radial velocities (RVs) with a 
precision of a few hundred m s“^, on the basis of previous 
experience with the instrument. 

A total of 2771 stars covering GMD positions from the 
upper red-giant branch to about one magnitude fainter than 
the main-sequence turnoff (TO) luminosity were observed 
between 2003 and 2009, including 306 new spectra obtained 
in 2009 and targeting MS stars already observed in the pre- 
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Table 1. Spectroscopic Observation Statistics 


No. Observations 

No. Stars 

No. Total 

1 

83 

83 

2 

1722 

3444 

3 

501 

1503 

4 

321 

1284 

5 

72 

360 

6 

20 

120 

7 

1 

7 

8 

27 

216 

9 

10 

90 

10 

12 

120 

11 

1 

11 

12 

1 

12 

Total 

2771 

7250 


vious epochs. Determination of the M4 velocity dispersion 
and binary star fraction were the prime motivators for ob¬ 
taining these data. Therefore nearly all stars were observed 
at least twice, and three or more spectra were obtained for 
nearly 40% of the sample. A total of 7250 individnal spectra 
were available for our study, as we summarize in Table 1. 

Default basic reductions (bias and dark correction, 
flat fielding) have been performed with the ESO pipeline. 
Spectral extraction has been performed using the opti¬ 
mal spectral extraction algorithm (Horne 1986) and the 
empirical point-spread function concept from photometry 
(Anderson & King 2000; Anderson et al. 2006) adapted to 
spectroscopy (LM14). The extracted spectra are fully con¬ 
sistent with the ones delivered by the ESO pipeline using 
either the HORNE option or the standard SUM, as a con¬ 
sequence of the fact that scattered light and fiber cross¬ 
talk are negligible in this specific setting of the MEDUSA 
mode. Wavelength calibration is performed by taking a hol¬ 
low cathode Thorium-Argon lamp (hereafter simply Th-Ar 
lamp) spectrum with each of the fibres in the MEDUSA 
plate with the selected instrument setup. These standard 
calibration frames are usually taken during the day, several 
hours before observations. During a science exposure five fi¬ 
bres called simultaneous calibration fibres (SimCals) can be 
dedicated to take additional Th-Ar spectra, in order to cor¬ 
rect for changes in the wavelength dispersion solution caused 
by variations in air pressure and temperature occurred in¬ 
side the spectrograph between the day-time calibration and 
night observations. S09 used the ESO Giraffe standard re¬ 
duction pipeline (Blecha et al. 2000) and the Ancillary Data 
Analysis Software (Royer et al. 2002) to reduce the M 4 data 
from raw CCD exposures to wavelength-calibrated ID spec¬ 
tra, apply the drift correction in the wavelength scale and 
determine the radial velocities of the observed stars. 

A systematic offset of 150 m s“^ in the radial velocity 
zero point was observed by in the five simultaneous calibra¬ 
tion fibres between 2003 and 2006 data (see §4.1 of S09 for 
a detailed discussion). Changes in environmental conditions 
inside the spectrograph (i. e. air pressure and temperature), 
as well as switching to different instruments modes, could 
cause an overall shift in wavelength (and hence RV) in the 
interval of time between observations and calibration expo¬ 
sures and between different nights. These should be auto¬ 


matically taken into account by the ESO pipeline using the 
simultaneous calibration fibres. After removing this source 
of systematic error, the only expected variation in RVs of 
calibration spectra taken in different epochs should have a 
random distribution, given by the precision in RV of the 
spectrograph and independent with time. 

The dispersion inside a single epoch is several times 
smaller then the offset (Fig. 6 of S09), and this seems to im¬ 
ply a problem in the original wavelength dispersion solution. 
Therefore we decided to revisit the basic reductions of the 
raw data sets, to see if the velocity precision could be signifi¬ 
cantly improved. In the next sections our different approach 
to wavelength calibration and drift correction is described. 
To mark the difference with the ESO pipeline, we will refer 
to our approach as Direct Calibration. 


3 GIRAFFE WAVELENGTH CALIBRATION 

ESO provides a pipeline for GIRAFFE as a part of the 
VLT Data Flow System (DFS). A detailed description of the 
pipeline can be found at the ESO website^; here we provide 
a brief description of the wavelength calibration subroutine 
to highlight the differences with our Direct Calibration tech¬ 
nique. 

3.1 ESO Pipeline calibration 

Calibration frames used for dispersion solutions are obtained 
by illuminating the entrance of the fibres with a hollow cath¬ 
ode Th-Ar lamp using the same setup as for the stellar spec¬ 
tra. A complete GIRAFFE dispersion solution consists of an 
optical model of the fibre spectra onto the CCD detector, a 
polynomial fit of the optical model residuals, and the cor¬ 
rection of individual fibre offsets. 

The solution is determined iteratively, i. e. at each itera¬ 
tion the previous model is updated until the measured emis¬ 
sion lines have a negligible radial velocity with respect to the 
tabulated values. A previous dispersion solution is useful for 
rapid computations, but, if necessary, a new dispersion solu¬ 
tion can be computed from scratch using the average optical 
model that accompanies the pipeline software. 

In detail, for each wavelength calibration procedure, a 
set of unsaturated Th-Ar emission lines is selected accord¬ 
ing to the instrumental setup of the images. The position 
of these lines on the CCD is determined using an analyti¬ 
cal optical model. The residuals between these positions and 
the measured ones are then modeled with a 2D Chebyshev 
polynomial. The global wavelength solution is then the sum 
of the optical model and of the polynomial residuals. See 
Royer et al. (2002) for a detailed description of the ESO 
pipeline, and Bristow et al. (2010) on the use physical mod¬ 
els for instrument calibration. 

3.2 Direct Calibration 

A more classical, direct approach has been followed in our 
work. For each individual Th-Ar spectrum, a selected list 
of emission lines with known wavelengths are measured and 

^ http://www.eso.org/sci/software/pipelines/ 
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then a piecewise polynomial (represented as Basis spline, 
de Boor 1977) is fitted to determine the wavelength disper¬ 
sion solution. 

The algorithm presented here is designed to antomati- 
cally perform those steps that in other tools snch IRAF^ re¬ 
quire human interaction. Although our reference wavelength 
range and resolution are the ones defined by the HR9 set¬ 
ting of GIRAFFE, we want to keep our algorithm as gen¬ 
eral as possible for rapid adaptation to different wavelength 
GIRAFFE ranges or for spectra obtained with different in¬ 
struments. 


3.2.1 Line lists of Thorium and Argon lines 

Several published Th-Ar line lists are available in which lines 
have been selected according to their stability and strength 
from high resolution spectra; see for example Murphy et al. 
(2007); Lovis & Pepe (2007). However, when dealing with 
medium resolution spectroscopy, blending with close or faint 
lines becomes important. The centers of lines of interest can 
be shifted several hundredths of Angstroms, thus leading 
to a less precise determination of the wavelength dispersion 
even with excellent data and a very accurate line list. 

To create a line list that takes into account line blend¬ 
ing, we proceeded in the following way. We have selected 
the Th-Ar reference spectrum obtained with ESO’s 3.6m 
telescope Coude Echelle Spectrometer (GES, Enard 1982) 
as our reference among several Th-Ar spectrum available in 
the literature after checking the similarity with GIRAFFE 
Th-Ar spectra (emission features may change depending on 
the manufacturer of the lamp). This spectrum has been ob¬ 
tained with a resolving power of i? = 200,000, which is 
almost nine times the resolving power of GIRAFFE in the 
HR9 setting. We degraded this spectrum to match the reso¬ 
lution of our spectra; we will refer to this smoothed spectrum 
as the reference one. All the calibrations will be relative to 
the reference system given by this spectrum. 

Lines have been automatically identified using the first 
derivative of the reference spectrum. Fig. 1 shows an exam¬ 
ple of the line identification process, with the smoothed ref¬ 
erence spectrum in the upper panel (in arbitrary units /) and 
its first derivative in the lower one (in arbitrary flux units 
per Angstrom /A“^). To separate strong lines from small 
features and weak emission lines, we have selected only lines 
with null derivative at the expected central wavelength and 
inflection points on the side that have a minimum absolute 
value of 100 /A“^. 

In our wavelength range 178 lines have been identified. 
The positions of these lines have been measured on the refer¬ 
ence spectrum, and their values define our rest-frame system 
of reference for the determination of the wavelength disper¬ 
sion solution of GIRAFFE calibration spectra. 


^ IRAF is distributed by the National Optical Astronomy Ob¬ 
servatory, which is operated by the Association of Universities 
for Research in Astronomy (AURA) under cooperative agreement 
with the National Science Foundation. 




Figure 1. An example of the line identification on the Thorium- 
Argon reference spectrum, with the smoothed spectrum in the 
upper panel and its first derivative in the lower one. Very weak 
lines have been excluded in our analysis if the inflection points of 
their first derivative are within the two blue dashed lines. Selected 
lines are denoted by a red vertical line. 

3.3 First approximation of the wavelength 
dispersion solution 

A first approximation to the dispersion solution is needed to 
allow the algorithm to determine the solution. 

We start by selecting the ten lines with the largest fluxes 
from the reference Th-Ar spectrum in the wavelength range 
of interest. The number of lines has been chosen in a way 
that the included lines form a subset with mean flux that 
is at least a few times the flux of the brightest lines not 
included in the subset. This ensures that the selection in the 
observed spectra will not be affected by differences in the 
line flux ratios that may be present when comparing spectra 
obtained with lamps made by different manufacturers. A 
very rough pixel position for each line must be provided also: 
a single value is good enough for all the fibres, despite their 
shift in the cross-dispersion direction. The wavelength values 
are used as initial guesses for a Gaussian fit, so a precision of 
one half of the full-width-at-half-maximum (fwhm) is good 
enough. This selection can be easily performed by visual 
inspection directly at the reference spectrum, and must be 
done only once for a given wavelength range or instrument 
setup. This is the only human interaction required by the 
algorithm. 

A Gaussian fit on these pre-selected lines is performed 
on both the pixel space for the observed spectrum and in the 
wavelength space for the reference one. The values are sorted 
in increasing order, so that the values in the pixel space have 
now an associated wavelength. The only precaution to be 
taken at this step is to check that no strong lines are present 
close to the border of the wavelength range (e. g. 10 A from 
each limit for our instrument setting), otherwise an error 
in the association of lines with pixels can occur. In fact, as 
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a result of the optical design of GIRAFFE, the wavelength 
range of the spectra varies slowly across the CCD, with a 
maximum shift of about 10 A for the fibres close to the 
center of the sensor respect to the first one. 

Before proceeding to the first approximate determina¬ 
tion of a wavelength solution, we need an additional step 
to ensure that cosmic rays or uncorrected hot pixels have 
not been identified as spectral lines. By inspection of our 
data frames, we noticed that real spectral lines have a total 
shift of around 100 pixels in the cross-dispersion direction 
(i. e. across the fibres). The maximum shift with respect to 
the first fibre around the center of the frame is 50 pixels, so 
we set this value as a threshold for checking proposed line 
identifications in the various fibres. This steady shift in line 
positions in different fibres is a consequence of the optical 
design of the instrument and it varies smoothly across the 
sensor, so for each wavelength value we perform a low-degree 
polynomial interpolation of the pixel values versus their fibre 
number. The fit is performed a second time after the removal 
of outlier pixel values. The approximate wavelength solution 
is finally obtained interpolating the wavelength values versus 
the corresponding pixel values with a low-order polynomial 
(degree=3) 

3.4 Wavelength dispersion solution for a single 
frame 

Now that an approximate solution is available, the determi¬ 
nation of the final dispersion solution is an iterative process 
where the resulting precision is increased with each iteration. 
These are the steps that we followed: 

• at each iteration more faint lines are measured in the 
calibration frame and associated with lines in the reference 
spectrum. 

• Lines with fwhm higher then a given threshold, nomi¬ 
nally twice the instrumental fwhm at that wavelength, are 
excluded from the list; 

• The dispersion solution is computed again, and the or¬ 
der of the polynomial can be increased according to the num¬ 
ber of added lines. 

• Lines that have an interpolated wavelength that differs 
from the reference value over a given threshold are rejected 
and not included in next iterations; the threshold is initially 
very large and it is rescaled at each iteration. 

In the last two iterations, the lines are weighted ac¬ 
cording to their brightness and their width in both the cal¬ 
ibration spectra and the reference one. After several tests 
Equation 1 was found to provide the best combination of 
these parameters. 



where Wi is the weight associated to the line i, acai,i and 
Uref.i are the widths of the line in the calibration and refer¬ 
ence Th-Ar spectrum respectively, and At is the measured 
amplitude of the line. An example of the weight association 
as a function of wavelength is given in the upper panel of 
Fig. 2 for a randomly selected calibration frame. 

In the last two iterations, the polynomial function is 
substituted with a piecewise polynomial (in the form of Basis 
splines or B-splines). To fit the curve to our data we made 



5150 5200 5250 5300 5350 

Wavelength [A] 


25 50 75 100 125 

Fiber number 

Figure 2. Two quantities associated with pixel-wavelength cal¬ 
ibrations are plotted as functions of wavelength. In the upper 
panel, the weight associated to each line is shown. In the lower 
panel, the difference between the wavelength determined using 
the dispersion solution for the measured center of the line, and 
the reference wavelength from CES spectra is plotted. The stan¬ 
dard deviation of this difference is reported in the figure. Points 
from different fibres are color-coded accordingly. 

use of the EFC routine included in the SLATEC Common 
Mathematical Library®, a collection of mathematical and 
statistical routines written in Fortran 77. This procedure is 
repeated individually for each fibre. 

In the lower panel of Fig. 2 the difference between the 
wavelength obtained with the dispersion solution (from the 
position in the pixel space) and the reference wavelength 
from the CES spectrum is shown as a function of wavelength. 
The standard deviation of this difference, determined using 
all the lines measured in a calibration frame, is on average 
cr = 0.0033 A, which corresponds to ~ 1/15**' of a pixel. 

3.5 Overall calibration 

After the wavelength calibration of a frame, the lines used 
for the dispersion solution are stored. Our dataset span a 
large temporal range, from almost the beginning of the sci¬ 
entific observations of GIRAFFE in 2003 until 2009. The 
information on emission Th-Ar lines from all these frames 
is collected in order to identify a stable set of lines and im¬ 
prove the dispersion solution. Only the lines that have been 
successfully measured on 90% of the fibres are retained. The 
only exception is for those lines that are closer than 10 A to 
the edge of the spectrum, for the reasons explained in §3.3. 
The weight associated to each line is the mean of the weights 
calculated in each calibration fibre using Equation 1. A to¬ 
tal of 114 lines has been selected (Table 2). The dispersion 

® The official repository of Version 4.1 ( July 1993) is 
http://www.netlib.org/slatec/ 
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Table 2. Lines and Weights used for the Dispersion Solutions 


n 

Wavelength [A] 

Weight 

n 

Wavelength [A] 

Weight 


1 

5134.7457 

8.776 

58 

5239.5502 

8.574 

2 

5136.1202 

26.557 

59 

5241.8563 

30.608 

3 

5137.4775 

15.334 

60 

5243.7542 

24.384 

4 

5140.7642 

13.070 

61 

5247.6542 

9.685 

5 

5141.7838 

5.382 

62 

5250.8742 

18.822 

6 

5143.9137 

8.399 

63 

5252.8043 

11.150 

7 

5145.3047 

6.584 

64 

5254.2542 

12.568 

8 

5148.2122 

19.641 

65 

5258.3603 

4.827 

9 

5149.2102 

21.811 

66 

5260.1046 

9.866 

10 

5151.6104 

5.564 

67 

5261.4763 

15.428 

11 

5154.2434 

5.420 

68 

5262.6055 

23.083 

12 

5157.4512 

26.060 

69 

5264.8003 

19.409 

13 

5158.6042 

3.818 

70 

5266.7102 

6.025 

14 

5159.5911 

19.307 

71 

5269.7842 

19.385 

15 

5160.7159 

9.8942 

72 

5272.9343 

14.107 

16 

5161.5402 

4.598 

73 

5274.1183 

9.093 

17 

5162.2882 

6.315 

74 

5277.4967 

13.711 

18 

5163.4582 

6.360 

75 

5281.0698 

9.834 

19 

5164.4022 

23.736 

76 

5282.4010 

23.019 

20 

5165.7682 

12.135 

77 

5283.6902 

24.195 

21 

5166.6522 

23.825 

78 

5286.8895 

15.986 

22 

5168.9242 

18.375 

79 

5289.9048 

21.810 

23 

5170.2482 

23.561 

80 

5291.8202 

13.201 

24 

5170.2482 

25.057 

81 

5294.3981 

11.068 

25 

5173.6782 

14.890 

82 

5295.0866 

17.948 

26 

5175.3262 

11.583 

83 

5296.2823 

9.744 

27 

5176.9642 

6.212 

84 

5297.7575 

8.974 

28 

5178.4602 

21.545 

85 

5300.5243 

11.599 

29 

5180.7198 

22.175 

86 

5301.4042 

13.648 

30 

5182.5242 

15.943 

87 

5303.4689 

18.513 

31 

5184.4562 

19.374 

88 

5305.6923 

16.192 

32 

5186.4122 

16.782 

89 

5307.4660 

14.933 

33 

5187.3402 

5.066 

90 

5309.6143 

20.676 

34 

5190.8742 

19.680 

91 

5310.2690 

16.136 

35 

5194.4582 

9.513 

92 

5312.0023 

5.200 

36 

5195.8142 

5.877 

93 

5315.2242 

19.243 

37 

5197.2362 

21.669 

94 

5317.4943 

10.140 

38 

5199.1702 

4.668 
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21.339 
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21.904 
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20.316 
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17.428 
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100 
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46 
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24.140 

47 
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48 
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6.626 
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19.221 

49 

5220.9031 

8.875 
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50 
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22.528 

51 
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14.010 
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20.157 

52 

5228.9985 

22.407 
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5347.9737 

12.361 

53 

5231.1600 

3.836 

110 

5349.0053 

19.482 

54 
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18.568 
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5351.1289 

12.436 

55 
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19.372 
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21.831 

56 
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28.873 

113 

5353.0208 

22.987 

57 

5238.8142 

8.417 

114 

5354.6055 

17.678 


solution is then computed again, using the same function of 
the last iteration in §3.4 but with the improved line list. 

Using the same set of lines and weights to compute 
the dispersion solution results in a more homogeneous cal¬ 
ibration of the entire dataset, while using the more stable 
lines reduce the scatter around the computed dispersion. In 
Fig. 3, the resulting overall calibration for the same frame in 
Fig. 2 is shown. A further reduction of the dispersion from 
(T = 0.0033 A to (j = 0.0026 A of the difference between 
measured and reference wavelength of individual lines is ob¬ 
tained. To have an idea of the goodness of the dispersion 
solution, this result can be compared with the average pixel 
size of ~0.05 A, or typically < 1/19^^ of the size of a pixel. 



25 50 75 100 125 

Fiber number 

Figure 3. As in Fig. 2, on the same calibration frame but using 
the refined line list. Note the reduction in the standard deviation 
compared to the previous figure. 

4 SIMCAL CALIBRATION 

A simultaneous calibration lamp provides a reference sys¬ 
tem obtained simultaneously with the observations, to take 
into account drifts of the spectrum in the wavelength space 
that may occur during the time between the Th-Ar lamp 
calibration frame and observations. 

GIRAFFE has five fibres available for Th-Ar SimCal, 
homogeneously distributed between the science fibres. In the 
ESO pipeline, drifts are corrected by cross-correlating the 
SimCals with the provided numerical mask; the resulting 
offsets are linearly interpolated across the CCD in the cross¬ 
dispersion direction to determine the wavelength offset to 
be applied for each fibre to the dispersion solution obtained 
using the day-time calibration frame. Two apparently rea¬ 
sonable assumptions are made here: the drift is constant in 
the velocity space across the wavelength range of the consid¬ 
ered order, and that the overall shift is linear in the spatial 
direction of the CCD. 

We first check if these two assumptions are correct. In 
order to do so, we need to compare two Th-Ar calibration 
frames with a measurable drift. We then opted to compare 
calibration frames from two different nights. Emission line 
positions for each fibre are measured in both frames and the 
difference in pixel space is taken. Conversion of differences 
from pixel space into radial velocity space ARVi for 

the ith line is given by: 

ARVi = Axf *5A(f) c/\i (2) 

with Xi given by the day-time wavelength dispersion solu¬ 
tion, 5Xi is the size of the pixel in the wavelength space at 
the position of the ith line. 

Fig. 4 and the first panel of Fig. 5 show the difference 
between the calibration frames of two typical nights. Differ¬ 
ent night combinations show similar results, with the only 
significant difference being the amount of overall radial ve- 
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Figure 4. Radial velocity shifts between the calibration frames 
of two typical nights, along the cross-dispersion direction. The 
position along the wavelength dispersion direction is color-coded. 
This plot shows that the common practice of assuming a constant 
value along one axis is not justified. 


locity shift. It appears clear that both ESO pipeline assump¬ 
tions here do not hold: the shift is not constant for a single 
fibre, and it is not a linear function of the cross-dispersion 
direction. Probably these differential shifts also happen dur¬ 
ing individual nights, even if in a less pronounced way. We 
conclude that correcting for a single value of the radial ve¬ 
locity will not provide a good correction for drifts, specially 
when temperature and air pressure inside the spectrograph 
are not actively controlled, as in our case. 

In our new analysis, a 2D polynomial fit (Chebyshev 
polynomials of the first kind) of degree (2, 2) has been chosen 
to model the shift in radial velocity as a function of pixel 
coordinates. However, during science observations only five 
fibres out of 132 are available for simultaneous calibration 
(fibres 1, 32, 63, 94 and 125). To check if the polynomial fit 
works as well with a reduced number of fibres, we performed 
two fits of the measured RV differences for two given nights 
(first panel of Fig. 5): the first one using all the available 
fibres (second panel of Fig. 5), and the second one using 
only the SimCals fibres (third panel of Fig. 5), in both cases 
using the same function to fit the measured differences. 

In the example shown in Fig. 6, the difference between 
the drift model using only SimCal fibres and measured line 
positions (left panel) has a mean of (Ahv) < 10 m s“^ and 
standard deviation a ~ 125 m s“^ after 5 — a clipping. The 
difference with the model obtained using all the fibres (right 
panel of Fig. 6) shows a systematic of (Arv) < 10 m s“^ 
and a (T ~ 20 m s“^ (after 5 — a clipping). Considering that 
we are using only 5 fibres out of 132 to redetermine the drift 
across the whole detector, this is a reasonable result. Finally, 
we note that polynomials with higher degrees resulted in fit 
too sensitive to outliers when using only the five SimCal 
fibres. 


Figure 5. First panel: Radial velocity shift of Th-Ar emission 
lines between the calibration frames of two different nights. Sec¬ 
ond panel: 2-D polynomial fit to the data of the first panel, using 
all the available fibres in the frame. Third panel: the same fit is 
performed using only the Simultaneous Calibration fibres. Radial 
velocity shifts are color-coded. The chosen model is able to repro¬ 
duce the global shape of the RV shifts while using the SimCals 
alone, with little influence from outliers. 



Figure 6. First panel: differences between the SimCals drift 
model and the observed radial velocity shift. Second panel: dif¬ 
ferences between the SimCals drift model and the same model 
using all the calibration fibres. The same calibration frames of 
Fig. 4 and 5 have been used. The scatter in RV of individual 
lines is clearly larger than the error we introduce in the global 
determination of the drift by using only the SimCals fibres. 
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5 RADIAL VELOCITY DETERMINATION 

In next section the three most popular techniques in radial 
velocity measurement are briefly described: Classical Cross- 
Correlation Function (CCF), Numerical CCF (nCCF), and 
Synthetic Template Matching. The steps for determination 
of radial velocities for each fibre within each exposure are 
then described in detail. 

Classical Cross Correlation function: This is the 
technique developed by Tonry & Davis (1979) and imple¬ 
mented in the RVSAO package (Kurtz & Mink 1998)"^. The 
observed spectrum is rebinned into a linear scale in logarith¬ 
mic wavelength, continuum normalized and then the Fourier 
power spectrum (FPS) is computed. The Fourier power spec¬ 
trum of a pre-selected template spectrum is computed too, 
and the two FPS are cross-correlated to determine the ra¬ 
dial velocity shift. In some modifications of the algorithm, 
the CCF is computed using the spectra directly, in the log- 
wavelength space. 

Numerical CCF: This is the technique introduced by 
Baranne et al. (1996) and improved by Pepe et al. (2002), 
also known as CORAVEL-type CCF. The observed spec¬ 
trum is cross-correlated in wavelength space with a numer¬ 
ical mask with non-negative values at the central position 
of the line, and zero otherwise. Each value of the mask, 
or weight, is determined using a synthetic spectrum, tak¬ 
ing into account the depth of the line and its width, in order 
to optimize the extraction of the doppler shift information. 
The cross correlation function is constructed by shifting the 
mask as a function of the Doppler velocity, and integrating 
the product of the mask with the observed spectrum. As a 
result of this definition, the value of the RV is given by the 
minimum in the cross-correlation function. To avoid confu¬ 
sion we will refer to this definition as nCCF This technique 
has proven its reliability in the discovery of several exoplan¬ 
ets, and achieves best results when a large number of lines 
is available (e. g., Pepe et al. 2004). 

Synthetic Template matching: A synthetic spec¬ 
trum to be used as rest-frame reference is generated using 
the atmospheric stellar parameters from photometry. The 
observed spectrum is shifted in the RV space and for every 
shift Vr the of the difference between the observed and 
the synthetic spectrum is computed. The value of Vr that 
gives the lowest is taken as the radial velocity of the star. 

The main disadvantage of Classical CCF technique is 
that the rebinning process (required for the FPS compu¬ 
tation) can introduce correlated noise in cases of very low 
SNR spectra. Since we want to avoid spectra rebinning in 
the first place this technique and further derivation have not 
been considered in the following analysis. 


5.1 Application of the numerical techniques 

We describe now our implementation of the Synthetic Tem¬ 
plate matching technique. To begin, the sky spectrum is re¬ 
moved from the observed one and the result is normalized 
to the continuum level of the stellar flux as later described 
in §5.2. Instead of a x^ function, we minimize the least- 
square function LSQ defined in Equation 3, where /* is the 




Figure 7. Comparison between the “CORAVEL-type” Numeri¬ 
cal CCF (upper panel) and Synthetic Template Matching (lower 
panel) on the spectrum of a SGB star with V = 16.4. The blue 
line represents Gaussian fits to these functions. 


observed stellar spectrum (after the steps described above), 
/syn is the synthetic one, and the sum is over the data 

points of the observed spectrum. Pixels affected by cosmic 
rays or CCD defects and spectral ranges contaminated by 
Th-Ar emission from close SimCals fibres are flagged with a 
zero value in mask, and unitary value otherwise. 


'^pixel 

LSQ{Vr)= 


/*(^) fsynji, Vr) 
1.2 -f fsyni^i, Ut") 


2 

mask{i). 


(3) 


For each determination of the LSQ function the synthe¬ 
sis is shifted by —Vr in the radial velocity space (in the oppo¬ 
site direction since Vr is the velocity associated to the star) 
and then rebinned into the wavelength scale of /*. This ap¬ 
proach avoids the more problematic rebinning of the noisier 
observed spectrum, and is the reason why the dependence on 
Vr has been included in the synthetic spectrum /syn(*, —Vr)- 
The denominator of Equation 3 differs significantly from 
the one of a standard x^ function. We did not include the de¬ 
pendence on the continuum in LSQ determination because 
this shape is changing with time and changes in instrumen¬ 
tation (e. g. the change of CCD in May 2008®); thus it could 
have introduced a bias depending on the epoch of observa¬ 
tions. We instead decided to divide by the synthetic spec¬ 
trum to give more weight to spectral lines, while a constant 
has been added to avoid an excessive weight for the deepest 
lines. The value of the additive constant has been found em¬ 
pirically by performing several tests on a subsample of stars 
of different spectral kind and with multiple observations. 


® See the FLAMES Giraffe Data Pro- 
A detailed description of the package is available at cessing and Quality Control website 

http://tdc-www.harvard.edu/iraf/rvsao/xcsao/xcsao.proc.html http://www.eso.org/observing/dfo/quality/index_giraffe.html 
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We have tested the Numerical CCF and the Synthetic 
Template Matching on several stars with spectra at various 
SNR to choose the technique that better suits our dataset. 
The resulting functions for a SGB star with V = 16.4 are 
shown in Fig. 7 in order to show the differences between 
the two techniques. For a given value of the radial velocity 
the Numerical CCF technique is using only the portions of 
the spectrum that are falling inside the positive values of 
the mask (namely <5—functions), while with the Synthetic 
Template Matching the LSQ function is determined over 
the entire spectrum. In both cases the width of the result¬ 
ing functions will be given by the sum in quadrature of the 
average width of the spectral lines in the observed spec¬ 
trum and the width of the lines in the mask/synthesis, hence 
with the Synthetic Match technique the final functions will 
be broader with respect to the Numerical derived one. A 
broader function will result in a greater uncertainty in the 
determination of its center, i. e. the Synthetic Match delivers 
a less precise radial velocity than the CORAVEL technique 
(if the same linelist is used in both cases). On the other 
hand, while this may be true in an ideal situation, when the 
limited wavelength range of the spectra allows to use only a 
few hundred of lines (with respect to the thousands of lines 
as in Pepe et al. 2004) as in our case, the CORAVEL tech¬ 
nique is too sensitive to noise, as exemplified in Fig. 7. We 
therefore decided to use the Synthetic Match to determine 
the RV of the whole set of stars. 

5.2 Continuum normalization and Sky flux 
determination 

Continuum normalization is one of the most daunting tasks 
in stellar spectroscopy. Generally speaking, there are sev¬ 
eral ways to deal with the continuum. When a spectrum 
has extended wavelength coverage as in the case of Echelle 
spectra, and the aim is to fit synthetic spectrum to the data, 
the continuum is obtained by fitting pre-selected feature-less 
spectral regions around the one to be analysed (as done by 
Kirby, Guhathakurta & Sneden 2008, but many other ex¬ 
amples are available in literature). In other cases as in the 
case of Equivalent Width measurements the interest is fo¬ 
cused on individual lines and it is possible to use those sur¬ 
rounded by at least a few Angstroms of local continuum (see 
for example the criteria used by Sousa et al. 2007 to create 
their linelist). 

Our dataset however presents several major complica¬ 
tions. The spectral range is very short, AA = 214 A, and 
very crowded, so there are no line-free regions available to 
set a continuum level reliably. In fact, the HR9 setting was 
chosen for the richness in this spectral window. The wings 
of the Magnesium-I triplet (5167.33 A, 5172.70 A, 5183.62 
A) affect nearly one quarter of the spectrum beginning at 
the low-wavelength edge of our spectra, so that the con- 
tinua must be extrapolated from the last three quarters of 
the spectra. Moreover, many lines are blended or very close, 
making the determination of local continua very difficult. 

To determine the continuum for each spectrum, we used 
a technique similar to the one introduced in LM14. Briefly, 
the continuum level and sky contribution are determined 
by using a spectral synthesis, obtained using the photomet¬ 
ric stellar parameters and average metallicity of the clus¬ 
ter, and the Solar Flux Atlas (Kurucz et al. 1984) as a 


template sky spectrum, both normalized to unity. Photo¬ 
metric atmospheric parameters for M4 have been derived 
using the photometric catalog kindly provided by Peter 
B. Stetson (D’Antona et al. 2009) and the calibrations from 
Ramirez & Melendez (2005) and Casagrande et al. (2010). 
The latest determinations for reddening and distance from 
Hendricks et al. (2012) have been adopted: E{B — V) = 
0.37 ± 0.01, Rv = 3.62 and Dciuster = 1.80 kpc. Synthe¬ 
ses are generated using the current version of the Local 
Thermodynamic Equilibrium (LTE) code MOOG (Sneden 
1973), and Kurucz (1992) atmosphere models calculated 
with the tr—enhanced New Opacity Function Distribution 
(AODFNEW, Castelli & Kurucz 2004). More details on the 
photometric parameter determination, atomic line parame¬ 
ters and synthesis calculation can be found in LM14. 

To apply this technique, we need to know the radial 
velocity of the star. At this stage high precision in RV is not 
required. A single pixel has an average size of 2.5kms“^, 
so an error of several hundreds m s“^ still yields a good 
fit of the synthesis to the observed spectra, given the fact 
that our goal at this stage is to determine the continuum 
level and not to derive accurate radial velocities or stellar 
parameters. The radial velocity is derived using a numerical 
CCF (Baranne et al. 1996) with line positions and weight 
derived from the same synthetic spectrum that will be used 
to determine the continuum, following the prescriptions in 
Pepe et al. (2002). 

i i r\i +S\/2 

nCCF{vr)^^nCCF\vr) = ^ " f{\yd\. 


As can be seen from Equation 4, the value of the nCCF 
at a given radial velocity point Vr is proportional to the 
sum of the stellar flux collected by each 5 function i in the 
numerical mask, rescaled for the weight a. In this equation, 
A(,^ is the wavelength of the line i shifted for the radial 
velocity of the nCCF and 5\ is the size of the hole in the 
numerical mask. 

The result of the numerical CCF allows us to deter¬ 
mine the star/sky flux ratio. The nCCF is computed for 
the smoothed synthetic spectrum (nCCFsynth) and the sky 
spectrum (nCCFsky) without any further correction and us¬ 
ing the same numerical mask. The cross-correlation function 
of the observed spectrum nCCFstar is then fitted with a 
linear combination of nCCFsynth and nCCFsky, being the 
linear coefficients of the fractions of stellar and sky fluxes. 

The upper panel of Fig. 8 shows the nCCFsynth (red 
line) and the nCCFsky (blue line) derived from normalized 
spectrum. The sky spectrum has a lower continuum and a 
deeper nCCF because of the higher metallicity and lower 
temperature respect to the observed star. In the lower panel, 
the two nCCFs are rescaled according to their estimated 
fluxes in order to show their relative influence, and their 
sum (green line) is compared with the observed nCCF (black 
points). Note here that this technique works best when the 
observed RV of the star (including Earth motion contribute) 
differs significantly from the RV of the sky, i. e. > 

20 kms“^. 

Both solar and synthetic spectra then are rebinned into 
the wavelength scale of the science spectrum (keeping con- 
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Figure 8. In the top panel, the nCCF from the normalized spec¬ 
trum of the synthesis (blue line) and the sky (red) are shown. In 
the bottom panel, nCCFsyn and nCCFgi^y are rescaled accord¬ 
ing to the measured flux to show the relative influence, and their 
sum (green) is compared with the observed nCCF (black). A star 
with Vmag = 16.4 Is Considered here as example. 

stant the flux for wavelength unit) and degraded to match 
the spectral resolution of the instrument. Finally the in¬ 
strumental effects (e. g.flat-fielding) are applied to make the 
comparison with the observed spectrum possible. A normal¬ 
ized linear combination of the two modified spectra, with 
the sky/fiux ratio as coefficients, is multiplied by a polyno¬ 
mial function: its coefficients are determined by fitting the 
observed spectrum. The partial continuum levels of star and 
sky are obtained by substituting the initial synthesis and 
the sky spectra respectively with unitary spectra and apply¬ 
ing all these same steps, with the exception of the fit of the 
polynomial function which coefficients are now known. The 
overall continuum level is the linear combination of the two 
partial continua. 

When sky subtraction and continuum normalization 
play a decisive role, e. g. in the determination of the stellar 
atmosphere parameters, more refined values for the star/sky 
ratio and continuum coefficients should be obtained as de¬ 
scribed in LM14, where a much larger emphasis was given 
to the quality of the sky subtraction and continuum normal¬ 
ization. In that case, the star/sky ratio was fitted simulta¬ 
neously with the stellar parameters of the synthetic spectra 
and the continuum placement, while here it is determined 
prior to the RV measurement. 

Atmospheric parameters and abundance measurements 
are particularly sensitive functions of the depth of the spec¬ 
tral lines, e. g. underestimating the continuum level would 
result in lower abundances, so refined values for the star/sky 
ratio and continuum coefficients are required. Compared to 
this, measuring the shift in wavelength of a spectrum is a 
simpler process. Thus while removing the sky contribution 
and removing instrument signatures can help, ultimately the 
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Figure 9. Upper panel: the contribution of sky (brown line) and 
star (red line), their combination (orange line) and the derived 
normalization level (yellow line) compared to the observed spec¬ 
trum (black points) for a program star with Vmag = 16.4. Middle 
panel: the same spectra zoomed in on the Magnesium triplet re¬ 
gion. Lower Panel: the residuals of the fit (observed spectrum 
minus green line) have a standard deviation of 21 e“, consistent 
with the SNR of the spectrum. 


achievable precision will depend on other factors such as the 
quality of the wavelength calibration and the SNR of the 
spectra. For this reason we have decided to stick with the 
simpler approach presented here. 

Fig. 9 shows the continuum determination for a Vmag = 
16.4 star, with photometrically derived Tes = 6000 K, 
log g = 4.0, and [Fe/H] = —1.15. The spectrum was taken 
on 2006 September 04, about five days before full Moon, 
at which time the angular separation between M4 and the 
moon was ~ 40°. Extraction of a pure stellar spectrum here 
is very difficult: the SNR is low, the sky contamination is 
large, the stellar lines are strong, the Mg I triplet in star 
and sky spectra severely depress the fluxes at the blue end 
of the spectral range, and as always the total spectral range 
is short. All these effects combine together to make it al¬ 
most impossible to find suitable line-free continuum zones 
to normalize the spectra. However, as we show in the lower 
panel of Fig. 9, our technique works well even in these diffi¬ 
cult cases. The observed spectrum is closely matched by the 
combination (green line) of the contribute of the sky (red 
line) and the stellar spectrum (blue line), with a sky/star 
ratio of 0.46. The knowledge of the individual contributes 
allows a reliable continuum determination even in problem¬ 
atic region of the spectra, such as around the Magnesium 
triplet. The residuals of the fit in this region (in the lower 
panel) have a standard deviation of 21 e~, which is very 
close to the value of the photon noise, ~ 17e“. 

For faint stars, even if observed only a few days after 
new moon, the sky can contribute up to half of the observed 
spectra of faint M4 targets. Otherwise the additional lines 
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Figure 10. The LSQ function for the same observations shown 
previously, before (red line) and after (black) sky subtraction. 
The blue line represents the gaussian fit to determine the RV of 
the star. 
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Figure 11. Distribution of the radial velocity dispersion (t„ as¬ 
sociated to each star before (red lines) and after (black lines) the 
SimCals correction. Vertical lines mark the median values of the 
two distributions. Three ranges in de-reddened V magnitude are 
considered. 


would lower the continuum level determined using only the 
stellar synthesis. Fortunately, our spectral range (5140A < 
A < 5360A) is free from atmospheric absorption lines, so no 
correction is required for this potential issue. 

Fig. 10 shows the LSQ function for this observation 
before and after the removal of the sky spectrum. The spec¬ 
trum has been normalized using the continuum associated 
to the combined sky-|-synthesis in the first case, and the 
continuum associated to the synthesis alone in the second 
case after sky removal. The slope introduced by the contri¬ 
bution of the sky, which inevitably alters the shape of the 
curve associated to the star, is essentially gone after sky- 
subtraction from the observed spectrum. Finally, confidence 
in this approach for removing sky contamination and con¬ 
tinuum normalization is increased by noting that stars with 
multi-epoch observations show lower internal dispersion in 
their RVs than can be obtained with other methods. 


6 RESULTS 

After sky subtraction and normalization, the radial velocity 
is finally determined with a gaussian ht of the resulting LSQ 
function, determined as described in 5. 

Radial velocities shifts due to the intrinsic velocity of 
the star and the Earth motion are then applied to the stellar 
synthetic spectrum. Heliocentric and Barycentric corrections 
are determined using the FORTRAN subroutine by Stumpff 
(1980). 

The error in the radial velocity associated with a single 
exposure is a function of the spectral type of the star, the 
characteristics of the instrument, the radial velocity deter¬ 
mination technique and the SNR of the spectrum. Instead 
of trying to determine a formal velocity error of a single ex¬ 


posure, we determine an empirical error using the several 
radial velocity measurements available for each star. 

The weighted mean velocity of a star Vr is determined 
using Equation 5, where r is the index of the star, j is the 
index referring to individual exposures and n is the total 
number of observations available for the given star. 


E j — 1 '^T,j 


(5) 


We use as weight Wr,j the inverse of the minimum of the 
LSQ function (Fig. 10) after normalizing it for the number 
of pixels used. 

Following S09, the unbiased estimator for the popula¬ 
tion variance is determined with 


Cvr) = 


E li 




j — UJr j Dj' 


( 6 ) 


It is interesting to compare the distribution of the radial 
velocity dispersion associated to each star before and af¬ 
ter the SimCals correction for three different intervals in 
de-reddened V magnitude (Fig. 11). Going from the lower 
panel to the upper one, we can see how the correction for in¬ 
strumental drifts in the radial velocity becomes increasingly 
more important for brighter stars, i. e. with higher SNR. This 
improvement is highlighted by the difference in the median 
values of the two distributions (vertical lines in the Figure). 
Although radial velocity measurements for faint stars are 
mostly dominated by photon noise, a small improvement in 
the cr„,. is still present. 

The average scatter (u^,,)*’*" for a magnitude bin is com¬ 
puted as the 68.27'^*' percentile of (t„,. distribution. This value 
is computed for several magnitude bins and the resulting val¬ 
ues are interpolated to determine the estimated error in the 
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Figure 12. The unbiased weighted radial velocity dispersion de¬ 
termined for stars with multi-epoch measurements. The estimated 
error on the radial velocity as a function of magnitude {<yv,.){V) 
from our analysis (red thick line), is compared with the one ob¬ 
tained by Sommariva et al. (2009) (blue line) using the same 
dataset. The histograms of the distribution for different ranges 
of magnitude are shown on the right panels, where the vertical 
lines follow the same convention defined in the left panel. Stars 
that satisfy the condition > 3(crv^) are considered binary 
candidates. 


radial velocity as a function of magnitude (cr„,.)(y) (the sub¬ 
script is omitted for the sake of brevity in the next sections). 
In the left-hand panel of Fig. 12 we show the computed (cr„,.) 
values for our M 4 stars that have multi-epoch observations. 
In addition, we show the comparable trend that was derived 
by S09 (red thick line in Fig. 10 of their paper, blue line in 
Fig. 12). In the right-hand panel of this hgure we show his¬ 
tograms of in hve magnitude bins. To ensure a fair com¬ 
parison, we did not include the data acquired in 2009 in this 
analysis (see Section 2). At high SNR (i. e. Vb < 13) the er¬ 
ror is dominated by the instrumental stability, the goodness 
of the SimCals correction and the precision of the RV deter¬ 
mination technique, for a total uncertainty of ~ 100 m s“^ 
(< 0.0018 A, or < 0.02 pixel). We regard this as the prac¬ 
tical precision limit in radial velocity that can be obtained 
from these GIRAFFE data. 

For brighter (RGB) M 4 stars the effect of the sky in RV 
determination is negligible even for observations obtained 
close to the full moon. For fainter stars with lower SNR data, 
the quality of the spectra becomes the main component in 
the precision of the RV determination. The combination of 
an improved wavelength calibration and careful determina¬ 
tions of radial velocities result in an increased precision of 
about a factor of two compared to the S09 results in the 
radial velocity determination of a single star. 

Following S09, we consider stars with > 3(cr„,.} to 
be M4 binary candidates. Using their same dataset, i. e., 
without data taken in 2009, we find 22 candidates out of 
454 targets with Vb < 14.9, and 55 candidates out of 2068 


targets with Vb > 14.9, for a total of 77 binary candidates. 
A total of 20 more candidates are found compared to the 
previous analyses by S09, using the same dataset and the 
same statistical approach. When the data from 2009 are in¬ 
cluded in the analysis 10 additional binary candidates with 
Vb > 14.9 are found, with the total number of candidates 
growing to 87. Results are summarized in Table 3, while 
individual RVs are reported in table 4. 

Target stars have been carefully selected using multi¬ 
band photometry and proper motions, drastically reducing 
the chance of contamination by other stars. If present, con¬ 
taminated spectra would have caused anomalous values in 
the FWHM or depth of the fitting function, which however 
have not been observed. 

We did not attempt to determine the cluster binary 
fraction because the limited number of available epochs re¬ 
quires detailed completeness simulations, which are beyond 
the scope of our work, but it is likely that the true fraction in 
the outskirts of the cluster is substantially less than 0.1, in 
agreement with Milone et al. (2012) and Nascimbeni et al. 
(2014, Paper III). A future paper of the series will deal with 
the fraction and spatial distribution of binaries in M4 with 
a joint analysis of photometry, RV and proper motions. 


6.1 Mean radial velocity and velocity dispersion 

The average velocity of M 4 has been computed for four mag¬ 
nitude bins. Binary candidates are excluded from the sam¬ 
ple, while we retain stars with only one RV measurement 
by assigning them (t„. = (cr„,.) evaluated at their magni¬ 
tude. Outliers in the CMD and obvious non-members (i. e. 
those stars with RV lower than 55 kms“^ or higher than 
87 kms“^) are excluded as well, leading to a total of 2543 
stars out of the 2771 of the initial sample. The weighted ra¬ 
dial velocity mean for the cluster, the associated error, the 
dispersion of the distribution, and the number of stars for 
each magnitude bin are reported in Table 5 and displayed 
in Fig. 13. Errors on RV dispersions have been calculated as 
described in § 6.2 There is no significant trend in the cluster 
mean radial velocity with magnitude. 

Our measured mean radial velocity for M4 is (vc) = 
71.08 ± 0.08 kms“^ with a cluster dispersion of ^vc = 
3.97 ± 0.05 kms“^. This value is in agreement with previ¬ 
ous studies. P95 reported a radial velocity mean of {v) = 
70.9 ± 0.6 kms“^ using 200 giant stars. Cote & Fischer 
(1996) derived a value of 70.3T0.7 kms“^ from the analysis 
33 turn-off dwarf stars. Our analysis marginally agrees with 
the one performed by S09, which determined a mean radial 
velocity of 70.27 ± 0.19 kms“^ using our same dataset. The 
main differences with that study, as described in previous 
sections, are in the wavelength calibration and in their use 
of a numerical mask. Unfortunately no radial velocity stan¬ 
dards are present in our dataset, so for the moment it is not 
possible to assess the origin of this small offset. 

Before discussing the RV dispersion as a function of the 
distance from cluster center, note from the two upper-right 
panels in Fig. 13 that the RV distribution for RGB and SGB 
stars deviates from a simple gaussian distribution. To verify 
that these deviations do not affect the determination of the 
radial dispersion profile of the cluster, in Fig 14 we have 
plotted the RV distribution of RGB stars (Vb < 14.75) for 
seven annuli of radial distance and for the full dataset (in the 
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Table 3. Average radial velocities of the 2791 stars of the sample. Only a sample is given here, the full table will be available in the 
online version of the manuscript, together with a table comprising the binary candidates only. 


a (J2000) 

S (J2000) 

ID 

V 

(B-V) 

Fo 

(B - V)o 

No. Obs. 

Vr 

km 

<Tvr- 

km 

{(Tvr) 

kms“^ 

(^v-r / i^Vr) 

kms~^ 

245.838250 

-26.513639 

29622 

16.329 

0.970 

14.974 

0.615 

2 

70.553 

0.115 

0.221 

0.521 

245.892417 

-26.504083 

40918 

14.505 

1.139 

13.146 

0.771 

2 

70.963 

0.023 

0.088 

0.266 

245.913167 

-26.483389 

46022 

15.016 

1.103 

13.649 

0.742 

3 

63.619 

0.084 

0.094 

0.891 

245.941333 

-26.560778 

51878 

12.935 

1.282 

11.570 

0.930 

10 

67.314 

0.266 

0.069 

3.834 

245.897083 

-26.537944 

42121 

12.944 

1.245 

11.566 

0.873 

10 

79.513 

6.666 

0.069 

96.145 


Table 4. Individual radial velocities of the 2791 stars of the sample. Only a sample is given here, the full table will be available in the 
online version of the manuscript. 


ID 

MJD 

j 

GIRAFFE File 

Fib 

BERV 

^obs 

r,j 

Vj-^j 

FWHM 

SNR 


34087 

52819.117744 

1 

GIRAF.2003-06-29T02:49:33.104 

2 

-13.951 

83.002 

69.051 

9.971 

7.25 

105.998 

34087 

53982.999641 

2 

GIRAF.2006-09-04T23:59:28.940 

22 

-29.533 

99.591 

70.058 

9.226 

12.08 

171.798 

30327 

52819.117744 

1 

GIRAF.2003-06-29T02:49:33.104 

3 

-13.957 

86.523 

72.566 

9.048 

8.61 

134.988 

30327 

53947.065775 

2 

GIRAF.2006-07-31T01:34:42.930 

5 

-25.373 

98.630 

73.257 

8.734 

16.38 

185.676 

30327 

53981.986284 

3 

GIRAF.2006-09-03T23:40:14.963 

5 

-29.544 

102.366 

72.822 

9.420 

12.66 

161.065 


Table 5. Radial velocity means 


N 


Vq range 


kms~^ 


km s“^ 


f^Vc 

km 


km s“ 


113 

168 

911 

1351 

totals: 

2543 


10.10 ^ Vo < 12.8 71.19 0.33 3.56 0.24 

12.80 ^ Vo < 14.15 70.97 0.33 4.33 0.24 

14.15 C Vo < 15.50 70.86 0.13 3.97 0.09 

15.50 ^ Vo < 16.85 71.26 0.11 3.94 0.08 

10.2 ^ Vo < 17.0 71.08 0.08 3.97 0.05 


last panel on the right). The gaussian fit obtained using the 
full RGB sample is shown for comparison in all the panels, 
after properly rescaling for the total number of stars in each 
sample. 

While the deviation from the normal distribution of 
each sample is not statistically significant, we have preferred 
to exclude the RGB and SGB stars from the following anal¬ 
ysis since in our judgement they cannot provide a reliable 
picture of the kinematical status of the cluster. It is also 
likely that evolved stars have different kinematics with re¬ 
spect to MS stars as a consequence of mass segregation, thus 
supporting our choice of analysing the MS sample alone. The 
faintness of MS stars fortunately is not a limit in the deter¬ 
mination of the kinematical properties of the cluster, thanks 
to the large size of the sample and our efforts in improving 
the RV precision. 


6.2 Radial velocity dispersion 

The dependence of the radial velocity dispersion with 
the distance from the cluster center is analysed by di¬ 
viding our sample into 10 radial, circular annuli on 


11 


12 


13 


14 


15 


16 



(B-V)„ 


V,, [km s“'] 


Figure 13. The weighted radial velocity mean of the cluster, the 
dispersion of the distribution and the number of stars for four 
magnitude bins, delimited with a grey line in the CMD. Candidate 
binary stars are denoted in red. There is no significant trend of 
the mean velocity of the cluster with magnitude. 


the sky projected (tangent) plane and centered on 
(q,< 5 )j 2 ooo.o=( 16 ‘' 23 “ 35 ! 22 ,- 26 ° 31 ' 32 '! 7 ) (Goldsbury et al. 
2010). The weighted velocity dispersion in each bin is then 
calculated, as displayed in Fig. 15. 

To take into account the broadening due to the obser¬ 
vational uncertainty, for each bin the 68.27**' percentile of 
the error distribution has been computed and the resulting 
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RV [km s"'] 


Figure 14. RV distribution of RGB stars (Vb < 14.75) for seven 
radial annuli (as defined in the label of each panel) and for the 
whole dataset (last panel on the right side). The distribution ob¬ 
tained by including all the RGB stars, rescaled for the number of 
stars inside each annulus, is shown for comparison (blue curve). 



Figure 15. The radial velocity dispersion is plotted as a function 
of distance from the cluster center (black circles). Thanks to the 
larger sample and the improved RV precision, our data show that 
the velocity dispersion is steadily decreasing as a function of ra¬ 
dius. This trend is suggested, but not statistically significant, by 
the data from Peterson, Rees & Cudworth (1995) (blue triangles) 


Table 6. Radial velocity dispersion as function of radial distance 


r 

arcmin 

C7 f 

arcmin 

f^Vc 

kms“^ 

Mint 

kms“^ 

^Mint 

kms“^ 

No. Stars 

0.5 

0.3 

4.50 

4.43 

0.52 

40 

1.5 

0.3 

4.53 

4.50 

0.24 

189 

2.5 

0.4 

4.21 

4.19 

0.19 

268 

3.5 

0.3 

4.14 

4.12 

0.18 

271 

4.5 

0.3 

3.97 

3.94 

0.19 

233 

5.5 

0.3 

3.79 

3.76 

0.19 

205 

6.5 

0.3 

3.74 

3.71 

0.20 

189 

7.5 

0.3 

3.36 

3.35 

0.20 

143 

8.5 

0.4 

3.19 

3.16 

0.21 

119 

10.5 

1.2 

3.48 

3.43 

0.24 

106 


value (~ 0.3 kms“^) quadratically subtracted from the ob¬ 
served cluster dispersion in each bin. The error in the 
mean has been quadratically subtracted as well. We use the 
symbol /tint to distinguish the resulting value from the ob¬ 
served cluster dispersion. The vertical error for each bin has 
been determined following Peterson & Latham (1986) for a 
consistent comparison with P95 measurements. The square 
in the uncertainty in /tfnt has been computed as the sum of 
sampling error = 2/Af/tfnt, where N the number of stars 
in each bin, plus the error due to the uncertainty in a single 
measurements ctJ = 4/(A'^ — Following the prop¬ 

agation of uncertainty, the error in /Tint is given by dividing 
the sum resulting from the previous step by 2/tint. Results 
are summarized in Table 6. 

Thanks to the larger sample and the improved RV preci¬ 
sion, our data show a steadily decreasing trend of the radial 
velocity dispersion as a function of radius. Our data are com¬ 
patible with the one of P95, although in their case the trend 
is not statistically significant. This trend is still present with 
the same amplitude under different assumed statistical re¬ 
strictions. These include: (a) considering only stars with at 
least three individual radial velocities; (b) applying a dif¬ 
ferent CT-clipping on RV errors, i.e. ^ 3{av^) instead of 
o'vr ^ 5(<^DT-)i (<^) selecting the stars in magnitude intervals. 
The radial distribution in P95 was determined using 182 
members with a radial velocity precision of ~ 1 kms“^ for 
individual stars, while our sample comprises a larger num¬ 
ber of stars (1763 MS members) and higher precision for a 
single star velocity determination (a few hundred m s“^). 

6.3 Rotation 

Cluster rotation has been determined by following the same 
approach of Cote et al. (1995). Briefly, the cluster is halved 
by the position angle PA, and the difference between the 
mean radial velocity of the two sides is computed, the er¬ 
ror associated to each point being the sum in quadrature 
of the error on the mean of each side. The rotation curve 
as a function of the PA is then obtained by repeating the 
measurement for different values of the PA. We adopted the 
same convention of Bellazzini et al. (2012; hereafter B12), 
i. e. the PA is increasing anti-clockwise from north (PA = 0°) 
to east (PA = 90°) in the plane of the sky, and the sine- 
function to fit the curve as AK- = Arot sin(PA-|-270 —PAq), 
where PAo corresponds to the PA of the rotation axis and 
Arot is the rotation amplitude projected on the plane of the 
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Figure 16. The difference between the mean velocities of the two 
halves of the cluster when divided by a line with varying PA, as 
a function of the PA itself. The best-fitting sine curve is shown. 
The analysis is performed on a sample including MS and SGB 
stars only (upper panel), and the same sample with only stars 
farther than 4' from the cluster center (lower panel). 


sky and averaged over the range of radius of the sample. 
To better investigate the properties of the cluster, we have 
performed the analysis twice, by including all MS and SGB 
stars (which should share the same dynamical properties), 
and then by considering only only stars farther than 4' from 
the cluster center, as previously done by P95. Binary candi¬ 
dates identified in Sec. 6 have been excluded from the anal¬ 
ysis. The best-fit values are A^ot = 0.51 ± 0.06 kms“^ and 
PAo = 127±7° for the whole sample, and Arot ~ 0.60 ±0.07 
kms“^ and PAo = 125 ± 7° for the distance-selected one, 
thus being perfectly consistent. (Fig. 16). These values are 
obtained using a Monte Carlo method. Each data point is 
perturbed of a random value extracted from a normal distri¬ 
bution with the associated error as standard rotation, and 
the rotation curve is fitted. This procedure is then iterated 
10.000 times to obtain a distribution of the two parameters 
and an estimate for their errors. 

A common test to assess the statistical significance of 
the rotation signal is to perform a two-side Kolmogorov- 
Smirnov test on the cumulative distribution of RVs of 
stars that he on opposite sides of the rotation axis (B12, 
Lardo et al. 2015). The presence of rotation induces a shift 
between the two distribution, see Fig. 17, hence the KS test 
quantify the probability to obtain the observed distribution 
under the null hypothesis that no rotation is present. We 
obtain Pks = 4.5% and Pks = 0.4% for the full sample and 
the distance-selected one respectively. The lower probability 
in the second case is expected, since we have removed from 
the sample the central region of the cluster, where the pro¬ 
jected rotational velocity becomes negligible. For the RGB 
stars alone we obtain Pks > 20 %, so we have not included 
the analysis in the text. Our results are consistent with the 



Figure 17. Comparison of the cumulative Vr distribution of stars 
on the left (dashed lines) and on the right (continuous lines) side 
of the rotation axis of the cluster. The cumulative distributions 
along with the probability of null-hypothesis of the Kolmogorov- 
Smirnov test Pks is shown for the whole MS-I-SGB sample (upper 
plot), and for MS+SGB stars farther than 4" from the cluster cen¬ 
ter. This result confirms the moderate rotation already observed 
in previous works. 

previous claim of moderate rotation for this cluster of P95. 
We found a rotation amplitude several times smaller than 
the value of Arot = 1.8T0.2 kms“^of Lane et al. (2010), al¬ 
though they quote half of this value as the rotation velocity. 


7 CONCLUSIONS 

We have presented a new approach for the determination of 
accurate wavelength dispersion solutions for not-stabilized 
instruments. We have developed an algorithm for improved 
corrections for radial velocity shifts using simultaneous cal¬ 
ibration fibres. Several precautions have to be taken when 
dealing with low SNR spectra covering a short wavelength 
range and strong absorption features, in particular during 
sky spectrum subtraction and continuum normalization. To¬ 
gether with a careful choice of the radial velocity technique 
that better suites the characteristics of our spectra, we have 
been able to significantly increase the precision of the ra¬ 
dial velocities with respect to existing tools. Algorithms have 
been optimized to require only a minimal user interaction. 

The code we have developed is lacking the necessary 
documentation to be user-friendly at present. While we do 
not exclude a general release in the next future, we have 
decided to not make it public yet. Instead, we provide a very 
detailed description of each step of our methodology with 
the hope that large consortium dedicated to data reduction, 
such as the Gaia-ESO or ESO itself, can implement and 
eventually improve it for other GIRAFFE settings as well 
other instruments. 

Our new radial velocity procedures have been tested on 
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a large dataset for globular cluster M4 with extant state- 
of-the-art radial velocity measurements in the literature, in 
order to assess the improvement in velocity precision. A to¬ 
tal of 7250 individual spectra for 2771 stars, gathered with 
the GIRAFFE spectrograph at VLT during a multi-epoch 
campaign, have been analysed, for a total of 2543 stars after 
removing non-members and candidate binaries. The derived 
cluster radial velocity is {v) = 71.08±0.08 kms“^, in agree¬ 
ment with previous measurements in literature but with an 
offset of An = 0.81 kms“^ with the previous value derived 
using the same dataset (probably due to the different radial 
velocity determination technique). We found a total of 87 
binary candidates, 22 giants and 65 sub-giants and dwarfs, 
which is a significative increase of the number found by pre¬ 
vious analyses. 

From our data we have measured a rotational amplitude 
of ~ 0.5km s~^, providing a statistical confirmation of the 
negligible rotation of this cluster. An average radial velocity 
dispersion of 4.5 kms“^ within 2' from the center of cluster 
and steadily decreasing outward has been found, in contrast 
with the nearly constant value of 3.5 kms“^ from P95. Our 
new determination is going in the same direction of recent 
simulations from Heggie (2014), where a N-body simulation 
of the cluster results in a RV dispersion steadily, and it will 
further improve the dynamical modeling of M4. Ultimately, 
the combination of these data with the unprecedented pre¬ 
cision of the proper-motions that will be delivered by the 
HST Large Program described in Paper I and following will 
provide a deeper insight into the three-dimensional dynam¬ 
ics of the cluster and a new, more accurate determination of 
its geometrical distance, among other things. 
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